Elimination of the linearization error in GW calculations based on the 
linearized augmented-plane-wave method 
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This paper investigates the influence of the basis set on the GW self-energy correction in the full- 
potential linearized augmented-plane-wave (LAPW) approach and similar linearized all-electron 
methods. A systematic improvement is achieved by including local orbitals that are deflned as 
second and higher energy derivatives of solutions to the radial scalar-relativistic Dirac equation and 
thus constitute a natural extension of the LAPW basis set. Within this approach linearization 
errors can be eliminated, and the basis set becomes complete. While the exchange contribution 
to the self-energy is little affected by the increased basis-set flexibility, the correlation contribution 
benefits from the better description of the unoccupied states, as do the quasiparticle energies. The 
resulting band gaps remain relatively unaffected, however; for Si we find an increase of 0.03 eV. 
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I. INTRODUCTION 

Electronic excitation energies may be obtained from 
the solution of the quasiparticle equation of many-body 
perturbation theory. This equation contains a nonlo- 
cal and frequency-dependent operator, the self-energy 
S'''^(r, r'; e), which, in principle, incorporates all elec- 
tronic exchange and correlation effects. As it cannot be 
treated exactly for real systems, practical implementa- 
tions typically use the GW approximation,^ which has 
become increasingly popular for electronic-structure cal- 
culations of excited states in recent years and yields band 
structures in good quantitative agreement with experi- 
mental spectroscopy for a wide range of materials. 

Due to its technical simplicity, the first implementa- 
tions were based on the pseudopotential plane-wave ap- 
proach. In spite of several approximations in the nu- 
merical treatment, which were necessary because of the 
lack of computer power in the 1980s, initial results were 
very promising. Hybertsen and Louie^ as well as Godby 
et al.^ showed that the calculated band gap of Si fell 
within a margin of about 0.1 eV from the experimen- 
tal value. Shortly afterwards the same authors reported 
band gaps for several other semiconducting materials 
that turned out to be equally accurate.'*'^ After these pi- 
oneering studies the GW approximation was applied to 
a variety of semiconductors, insulators, and metals with 
great success.^ 

So far, most codes still rely on the pseudopotential ap- 
proximation, which restricts the range of materials that 
can be examined. Transition-metal compounds and ox- 
ides, in particular, cannot be treated efficiently in this 
approach. Two early all-electron calculations using the 
GW approximation were done by Hamada et al7 for Si 
and by Aryasetiawan® for Ni, both within the linearized 
augmented-plane-wave (LAPW) method. However, only 



very recently were further full-potential implementations 
reported, based on the LAPW (Refs. 9 and 10), the lin- 
earized muf&n-tin orbital (LMTO) (Refs. 11 and 12), the 
projector-augmented- wave (PAW) (Refs. 13 and 14), and 
the Korringa-Kohn-Rostoker^^ method together with ap- 
plications to a larger variety of systems. 

Compared to earlier pseudopotential results, it was 
found, however, that LAPW, LMTO, and PAW calcula- 
tions appeared to yield systematically smaller band gaps 
for semiconductors and, in many cases, a worse agree- 
ment with experiment. ^'^^"^'^ Ku and Eguiluz^ hence ar- 
gued that, in contrast to all-electron methods, the pseu- 
dopotential calculations benefited from a fortuitious er- 
ror cancellation between the pseudopotential approxima- 
tion and the neglect of vertex corrections in the GW 
approximation. Although their calculations were subse- 
quently criticized for not being converged with respect to 
the number of bands, ^^'^^ other all-electron calculations 
showed a similar underestimation of the band gap.^^"^'* 
In a different attempt to make one step towards an all- 
electron treatment, Tiago et al}^ relaxed the pseudopo- 
tential approximation by constructing a pseudopotential 
only for the Is state of Si while treating the 2s and 2p 
states as valence. Surprisingly, the resulting band gaps 
did not deviate substantially from the previous pseudopo- 
tential results, which made the conjecture of Ku and 
Eguiluz doubtful. In order to resolve this conflict, it is 
imperative to carefully analyze and compare the numer- 
ical approximations made in the two approaches. 

Of course, deviations are expected for several reasons. 
First, the pseudized wave functions differ from their true 
counterparts and modify the matrix elements of the self- 
energy. Second, the core electrons are not included in 
the construction of the nonlocal self-energy if pseudopo- 
tentials are used, which may lead to errors in the core- 
valence exchange contribution. On the other hand, the 
single-particle wave functions (and the corresponding en- 
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ergies) of both pseudopotential and linearized all-electron 
approaches, such as LAPW or LMTO, become more and 
more inaccurate at higher energies. The LAPW basis set, 
on which we concentrate in the following, is defined by 
an expansion around fixed energy parameters that yield 
accurate wave functions only in their neighborhood, i.e., 
the valence band. Of course, this does not affect calcu- 
lations within density- functional theory^* (DFT), which 
only makes use of the occupied states. The GW self- 
energy, however, depends on the unoccupied states up to 
high energies through the Green function G as well as 
the screened Coulomb interaction W. The inappropriate 
description of these states might, therefore, cause errors 
in the self-energy correction. 

The construction of the pseudopotentials guarantees 
an accurate wave function and energy only for the ground 
state of a given angular momentum but leads to devia- 
tions for higher-lying states. While the deficiency in the 
pseudopotential approach is inherent in the pseudopo- 
tential construction, in the LAPW method it must be at- 
tributed to the inadequacy of the basis set for high-lying 
unoccupied states and can be overcome by increasing the 
basis-set flexibility. It is the purpose of this paper to elu- 
cidate the influence of the basis set on the GW results by 
systematically extending it towards basis-set complete- 
ness. We achieve this by adding local orbitals defined 
as second and higher energy derivatives of solutions of 
the radial scalar-relativistic Dirac equation. There are 
several alternative approaches: Bross and Fehrenbach^'' 
use spline functions to augment the basis set, Krasovskii 
et al.™ employ conventional local orbitals located in the 
conduction bands together with their energy derivatives. 
The advantage of the present approach is that no spe- 
cial consideration about the energy parameters is needed. 
Furthermore, with the order of the derivatives in each 
angular-momentum channel, it contains well defined con- 
vergence parameters that allow a systematic attainment 
of basis-set completeness. 

This paper is organized as follows. In Sec. II we provide 
a brief introduction to the GW and LAPW methods. In 
Sec. Ill we describe our extension of the LAPW basis set 
in detail. As an illustration, in Sec. IV the effect of the 
basis-set extension on the GW results for Si is discussed. 
We show that the extension of the basis set and the con- 
vergence with respect to the number of bands both re- 
duce the discrepancy with the pseudopotential and the 
experimental band gap. Unless stated otherwise, we use 
Hartree atomic units. 



II. METHODS 

A. GW approximation 

Within many-body perturbation theory, the quasipar- 
ticle wave functions ■(/'nk<7(r) and energies e„k<T are ob- 



tained from the solution of the quasiparticle equation 

(^-^V2 + F°'^*(r) + T/H(r)^ V^„k.(r) (1) 

+ / r'; e„k<T)V'nka(r') d^r' = e„k<7V'nk(7(r) , 

where F°'^*(r), V^{r), and (r, r'; e„kcr) are the ex- 
ternal potential created by the crystal field and other 
applied static fields, the Hartree potential, and the 
exchange-correlation self-energy, respectively. The quan- 
tum numbers n, k, and a signify the band, wave vec- 
tor, and spin. Like the majority of existing implemen- 
tations, our code exploits the formal similarity to the 
Kohn-Sham (KS) equation to obtain approximate ener- 
gies within first-order perturbation theory 



KS 
nka 



(2) 

where V^^{r) is the local exchange-correlation potential, 
¥'nkCT(r) the Kohn-Sham wave function, and the quasi- 
particle renormalization factor is given by 



Zrika — I 1 



de 



<1. (3) 



In the following we will always use the shorthand nota- 
tion (S-) =(^„k.|S-(eKSj|^„^^). 

We employ the GW approximation for the self-energy, 
symbolically written as Y,^'' = iGW, where G is the 
Kohn-Sham Green function 



G.(r,r';e) = ^ 



unocc 

+ E 

n,k 



(4) 



{5 is an infinitesimal positive number) and W the dynam- 
ically screened Coulomb interaction. The latter is calcu- 
lated from W = V + vPW, where v is the bare Coulomb 
potential and P the polarization function 



occ unocc 



n.k n'.k' 



^nka 



isy 



(5) 



x</'nk^(r)¥'«'kv(r)¥^;,k'<T(r')'/5nk^(r') 

in the random-phase approximation. In practice, the self- 
energy is decomposed into exchange and correlation con- 
tributions 



(6) 



The exchange contribution only depends on occupied 
states, whereas unoccupied states up to high energies, 

typically 100 — 200 eV above the Fermi energy, are needed 
for an accurate evaluation of the correlation contribution. 
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B. LAPW basis set 

The all-electron APW method^^ as well as the re- 
lated LAPW mcthod^^ rely on a decomposition of space 
into muffin-tin (MT) spheres, centered at the atomic nu- 
clei, and the interstitial region. The core-electron wave 
functions, which are (mostly) confined to the muffin-tin 
spheres, are directly obtained from a solution of the fully 
relativistic Dirac equation. Here only the spherical part 
of the effective potential is retained. For the valence elec- 
trons a basis set {(f>k+Gi^)} constructed. Its basis 
functions, the so-called augmented plane waves, are de- 
fined everywhere in space. The smoothness of the poten- 
tial in the interstitial region motivates the use of plane 
waves with |k + G| < -ftTnjaxj where -R'max is a convergence 
parameter. In the MT spheres, on the other hand, the po- 
tential is peaked at the nuclei and predominantly spher- 
ical. In the APW method, which is usually implemented 
in combination with a shape approximation of spherically 
symmetric potentials inside the MT spheres, one uses 
numerical solutions Ui^{ei„,r) of the radial Schrodinger 
equation 



with the Kohn-Sham Hamiltonian 



(0)/ 



1) 



2r2 



(7) 



(8) 



and the relevant effective potential V^'^ir). For simplic- 
ity we give the nonrelativistic equations here; the scalar- 
relativistic versions are shown in the Appendix. The 
atom index is suppressed throughout. Augmenting the 

interstitial plane waves with linear combinations of the 
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(e;o-,r) = u[J(e;o-,r)y;m(r), taking into account the 
continuity at the MT sphere boundaries, yields the ba- 
sis functions '/'k+cC'')- ^^"^ APW method the energy 
parameters are identical to the band energies enkcr- 
For a spherically symmetric MT potential this condition 
guarantees that the Kohn-Sham wave functions </5nko-(r) 
can be obtained exactly (subject to convergence with the 
Z cutoff and if max), but it leads to an expensive non- 
linear eigenvalue problem. In the LAPW method this 
difficulty is circumvented by fixing the e/o- at suitable en- 
ergies in the valence-band region, thereby avoiding the 
expensive self-consistency condition with respect to the 
band energies and making a linear solution of the Kohn- 
Sham eigenvalue equation possible. Apart from solutions 
U;°^(e;CT,'') of the radial Schrodinger equation (7), the 

,(0) 



dulJ{ei^,r)/de are also 



energy derivatives u\^{€i 
employed to increase the basis-set flexibility in the MT 
spheres. They are obtained from 



,(1) 



(0), 



de 



,(0) 



(0)/ 



(9b) 



The two sets of radial functions are combined with the 
interstitial plane waves in such a way that not only the 
basis functions but also their radial derivatives are con- 
tinuous at the MT sphere boundaries. This procedure 
yields the LAPW basis set for the valence electrons 



(r) 



i(k-|-G)T 



k+G,,(i^) 



if r ^ MT, 
(ei^,r) ifreMT, 



(10) 

where the coefficients a|^^ are uniquely determined by 

the matching conditions and is the unit-cell volume. 
The inclusion of the energy derivatives implies a linear 
approximation for the radial functions 



uia- (e„kcr, r) 



,(0) 



(e;o-,r) -I- (e„ko- - Q<t) 



(11) 

thus introducing a linearization error, which grows with 
the deviation of the band energy e„k<T from the param- 
eters ejo- and becomes especially relevant for high-lying 
states. In full-potential LAPW implementations the ra- 
dial functions are still derived from the spherical part of 
the MT potential. Therefore, the actual wave functions 
cannot be constructed from the wj^^^ (ej^- , r) alone, but as 
the nonspherical modulation is typically small, it is com- 
monly accepted that the energy derivatives add enough 
basis-set flexibility to describe the valence electrons in 
the full potential accurately.^^ 

In order to quantify the error incurred by the ap- 
proximation (11), we compare the Kohn-Sham eigenvalue 
spectra with and without the linearization of the radial 
functions. In the second case we determine the band en- 
ergies e„k(T iteratively by setting the energy parameters 
£(0- equal to the e„k(T until self-consistency is reached. For 
a spherically symmetric MT potential this procedure is 
equivalent to the APW method, because the contribution 
of the energy derivatives in our scheme vanishes at the 
self-consistency point eicr = f-nka, and the description of 
the wave functions is thus identical: linear combinations 
of 'U'[^^{eia,r) inside the MT spheres and plane waves 
in the interstitial region. It should be noted that al- 
though the APW basis functions exhibit a derivative dis- 
continuity at the MT sphere boundaries, the wave func- 
tions themselves are smooth. Furthermore, the increased 
basis-set flexibility achieved by the additional radial func- 
tions makes it possible to apply our method to the full 
crystal potential. In this case the energy derivatives yield 
a nonvanishing but small contribution that reflects the 
nonspherical modulation. 

In Fig. 1(a) we compare the resulting Kohn-Sham band 
structure with that from a standard full-potential LAPW 
calculation for the prototype semiconductor Si. Both 
band structures are evaluated within the local-density 
approximation for the same effective potential, which is 
obtained from the LAPW self-consistency loop. The cal- 
culation is carried out with the experimental lattice con- 
stant of 10.26 Bohr, a muffin-tin radius of 2.16 Bohr, 
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FIG. 1: Comparison of (a) full-potential LAPW and (b) psou- 
dopotontial band structures for silicon in the local-density 
approximation (dashed lines) with results from our nonlin- 
ear scheme in which the energy parameters are determined 
self-consistently (solid lines). In both cases deviations start 
to appear at around 20 eV above the Fermi energy (0 eV). 



^max = 4.0 Ha, an angular-momentum cutoff of 8, and 
512 k points in the full Brillouin zone. We use the 
FLEUR code.^'' In Fig. 1(b) we make a similar compar- 
ison with the pseudopotential plane-wave method, which 
is the most frequently used approach in the context of the 
GW approximation. For this calculation we use a stan- 
dard norm-conserving Hamann pseudopotentiaP^ and an 
energy cutoff of 9 Ha for the plane waves, the other pa- 
rameters are the same as in LAPW. In both cases devi- 
ations start to occur at around 20 eV above the Fermi 
energy, i.e., around the twentieth band. This is far below 
the number of bands normally included in a GW calcu- 
lation to ensure convergence of the self-energy. 

In the next section we discuss our extension of the 
LAPW basis set, which allows a systematic reduction 
and ultimate elimination of the linearization error. In 
the basis-set limit the basis becomes complete and all 
states, including high-lying conduction bands, are accu- 
rately described. 



III. THE LAPW BASIS-SET EXTENSION 

In order to enhance the basis-set flexibility, it is not 
sufficient to simply increase the plane- wave cutoff if max, 
as is the case in pseudopotential calculations, since this 
only improves the basis set in the interstitial region but 
not within the muffin-tin spheres. In fact, Krasovskii^*' 
has shown that a fairly small i^max is sufficient to provide 
enough flexibility in the interstitial region, while the MT 
part of the basis set quickly deteriorates the more the 
wave-function energies deviate from the parameters eia-. 



This inadequacy cannot be overcome with a higher i^max- 
In a straightforward extension of the LAPW approach 
one can include second energy derivatives in the Tay- 
lor expansion (11). In this case one would need a third 
matching condition, though, in order to uniquely define 
the augmented plane waves (10), e.g., the continuity of 
the second radial derivative at the MT sphere bound- 
aries. These more stringent conditions are known to lead 
to a less efficient basis set, i.e., one with a slower con- 
vergence with respect to i^Tmax, however. To circum- 
vent this problem we introduce the higher energy deriva- 
tives as additional basis functions in the form of local 
orbitals,^^ which vanish outside and on the MT sphere 
boundaries and, therefore, do not need to be matched 
to plane waves. In general, the radial part of a local 
orbital is constructed as a linear combination of three 
radial functions with the conditions of vanishing value 
and derivative at the MT sphere boundary as well as 
normalization. In our approach these three functions are 
u'il\ei^,r), u^ilHei^,r), and one ofu[''\ci^, r) with u>2. 
The latter is the vfh energy derivative, which can be ob- 
tained in analogy to Eq. (9) from 

(12a) 

^{n^i'J{^i^)WfJ{ei.))=Q (12b) 

(see the Appendix for the scalar-relativistic treatment). 

We now prove that the functions w^'^' (e, r) are (i) lin- 
early independent of each other and (ii) orthogonal to 
solutions u^^^{eo,r) of Eq. (7) with eo e that vanish 
outside and on the MT sphere boundary. For simplicity 
the indices I and a are omitted here. 

The first statement guarantees that each derivative in- 
creases the variational freedom and does not lead to an 
overcomplete basis set. For M(°)(e,r) and u^^\e,r) this 
follows trivially from condition (9b). Now assume that 
^(''^(e, r) for some v >2 can be written as a linear com- 
bination 

«('')(e,r) = ^c^«('')(e,r) (13) 

of the linearly independent fimctions iM'^ (e, r) with < 
/Lt < u—1. Then the apphcation of Eqs. (7), (9), and (12) 
leads to 

vru'^''-'^\e,r) = (h - ru'^''\e,r) 

= ^/xc^+iTO('')(e,r) , (14) 

which contradicts the assumption of linear independence 

of the u(^^(e,r) with < ^ < v—\. Therefore, the deriva- 
tives of orders 0, . . . , must all be linearly independent. 
The second statement guarantees that orbitals 

u^"^ (e, r) constructed from the valence basis are orthog- 
onal to the core states M*^°)(eo,r). It should be noted. 
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FIG. 2: Deviation of the Kohn-Sham eigenvalues from the 
results of our nonlinear reference calculation at the T point 

of Si as a function of energy for the basis sets (A), (B), (C), 
and (D) (for details see the text). The energy parameters are 
optimized for the valence bands and identical in all calcula- 
tions. The eigenvalues are given with respect to the Fermi 
level. The curves are smoothed using a Bzier algorithm for 
clarity. 



however, that shallow semicore states do not vanish suf- 
ficiently outside the MT spheres and should, therefore, 

be treated with (conventional) local orbitals at suitable 
energies. For the special case ly — the second statement 
follows from [h — €o)ru'^^\eo,r) = {h — e)ru^^\e,r) = 
with e 7^ eo and integration by parts, giving 



(e-eo) («(°)(eo)| u^'\e) 

/°)(eo,r) (hru^°\e,r)) 

- (W°)(eo,r-)) ru^°\e,r) 

^ (u(°)(e,/?)^z(°)'(eo,i?) -«(°Heo,i?)w'°"(e,i?) 




-f 

Jo 



dr 



(15) 



where R is the MT radius. For v > 1 this expression 

instead reads 



= ^ 
~ 2 



(e, R)u^°^'{eo, R) - u^"^ (eo, i?)u^"^'(e, R) 



(0)1 



,(-1) 







(16) 



where the last equality follows from induction. 

When applied to Si, the example shown in Fig. 1, 

already the inclusion of second-derivative local s, p, rf, 
and / orbitals yields agreement with the reference band 
structure from our nonlinear scheme, on the scale of the 
figure, up to 60 eV above the Fermi energy. The ad- 
dition of third derivatives and g functions pushes this 



limit up to around 80 eV. The improvement over conven- 
tional LAPW in the description of the unoccupied states 
is clearly seen in Fig. 2, which shows the deviation of 
the single-particle energies from the corresponding ref- 
erence values for (A) the conventional LAPW basis set, 
(B) with second-derivative local orbitals for I < 3, (C) 
with second-derivative local orbitals for / < 4, and (D) 
with second- and third-derivative local orbitals for Z < 4 
at the r point of Si as a function of the single-particle 
energies. The inclusion of just the second-derivative local 
orbitals for alU < 3 in (B) gives energies within 0.2 eV of 
the reference values up to 60 eV above the Fermi energy. 
Adding local orbitals with i = 4 in (C) and third deriva- 
tives in (D) further improves the agreement. The advan- 
tage over similar approaches to improve the flexibility of 
the basis set (e.g., Ref. 20) is that systematic convergence 
can be achieved without special assumptions about the 
energy parameters of the additional local orbitals. As we 
find the results to be relatively independent of their po- 
sition, one can simply use the LAPW energy parameters, 
which are located at the center of gravity of the valence 
band. The convergence towards basis-set completeness is 
then controlled by two simple parameters, the maximum 
quantum number I of the local orbitals and the order of 
derivatives u (together with the I cutoff and i^max of the 
augmented plane waves). 



IV. APPLICATION TO GW 

Once the wave functions are determined in the basis 
(10), the polarization function (5) and related quantities 
are represented in terms of a mixed basis designed for 
the expansion of products of eigenfunctions.^^ It consists 
of plane waves with a cutoff of 2.7 Bohr~^ in the inter- 
stitial region and products of two radial basis functions 
inside the MT spheres, where the first is related to an 
occupied and the second to an unoccupied state. For oc- 
cupied states contributions and higher energy deriva- 
tives can be neglected, since the parameter ei is chosen 
to be the center of gravity of the occupied Z-like density 
of states (eo = -7.9 eV, ei = -3.2 eV, ea = -3.4 eV, 
ei>3 = —4.8 eV relative to the top of the valence band). 
For unoccupied states the inclusion of with ly > 1 
as well as higher angular momenta can be important, 
because the corresponding wave-function coefficients be- 
come large for higher-lying states. In our calculations 
the mixed basis is constructed from uf^ with / < 2 for 

occupied as well as u'f^ with < Z < 6 and U;''^ with 
1< H<iyandl < L [(A) ly = 1, L = 3; (B) v = 2, L = 3; 
{C) p = 2, L = 4; (D) v = 3, L = 4] for unoccupied 
states. The omission of the other radial functions results 
in a negligible error not exceeding 0.1 meV in the cor- 
relation contribution (T,^) for all states considered here. 
The self-energy is evaluated with an angular-momentum 
cutoff of 6 and 216 k points in the full Brillouin zone. All 
other parameters are the same as in the underlying DFT 
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FIG. 3: Expectation value of the correlation contribution to 
the self-energy at the valence-band maximum of Si at F as 
a function of the number of bands. The inclusion of second- 
derivative local orbitals for Z < 3 (B) changes the asymptotic 
value by more than 0.05 eV with respect to the conventional 
LAPW basis set (A) . Further extensions of the basis give only 
minor corrections. 



TABLE 1: Expectation values of the correlation contribution 
{S'^) to the self-energy correction for the basis sets (A), (B), 
(C), and (D) (see text) together with the differences A(E'^) 
between basis sets (A) and (D). All values are in eV. 





(A) 


(B) 


(C) 


(D) 






0.646 


0.594 


0.585 


0.583 


-0.063 


Fl6c 


-4.183 


-4.224 


-4.230 


-4.233 


-0.050 


X41, 


1.824 


1.784 


1.778 


1.776 


-0.049 


Xlc 


-3.749 


-3.776 


-3.780 


-3.782 


-0.033 




1.139 


1.091 


1.083 


1.081 


-0.058 


Li. 


:!.<)li 


:!.<)(){) 


:!.<)()7 


:!.97J 


0.060 



TABLE II: Quasiparticle energies e for the basis sets (A), 

(B), (C), and (D) (sec text) together with the differences Ae 
between basis sets (A) and (D). All values are in eV. 





(A) 


(li) 




(U) 


Ac 




-0.646 


-0.685 


-0.692 


-0.694 


-0.048 


FlSc 


2.541 


2.514 


2.509 


2.509 


-0.032 


X4v 


-3.579 


-3.608 


-3.613 


-3.615 


-0.036 


Xlc 


0.521 


0.502 


0.499 


0.498 


-0.023 




-1.883 


-1.919 


-1.925 


-1.927 


-0.044 


Lie 


1.467 


1.432 


1.427 


1.424 


-0.043 



calculation. 

Changes in the; expectation values of the exchange term 
{H^) and the exchange-correlation potential {V^'^), which 
are both independent of the unoccupied states, are small 
and compensate each other; the main effect on the quasi- 
particle energies is due to the basis-set dependence of the 
correlation contribution to the self-energy. In Fig. 3 we 
show (H'^) for the valence-band maximum of Si at F as 
a function of the number of bands included in the Green 
fiuiction (4) and the polarization function (5) for the dif- 
ferent basis sets (see Sec. III). Up to the fiftieth band 
the curves are nearly identical but then begin to devi- 
ate due to the linearization error. As the denominator in 
(5) reduces the weight of higher-lying states while their 
linearization error grows at the same time, the largest 
increase of the difference between the curves is seen for 
intermediate numbers of bands between 60 and 150. The 
curves converge rapidly with respect to the basis-set size; 
already with the inclusion of second-derivative local or- 
bitals (B) convergence is reached to within 0.01 eV. 

The corresponding curves for other states look qual- 
itatively similar. Tables I and II give the values of 
(S*^) and the resulting quasiparticle energies e calculated 
with 245 bands for the different basis sets. As a ref- 
erence, the valence-band maximum in the underlying 
Kohn-Sham calculation is set to zero. All quasiparti- 
cle energies tend towards smaller values in the basis-set 
limit. Consequently, the effect on relative transition en- 



ergies is smaller but of the same order. In Fig. 4 we 
show the behavior of the indirect band gap as an ex- 
ample. The results are lowered by 0.02 eV if the k-point 
mesh is fully converged, and by another 0.02 eV if screen- 
ing due to the 2p electrons is included in the correlation 
self-energy. The convergence with respect to the num- 
ber of bands and the basis-set extension both increase 
the calculated band gap and thus narrow the distance to 
the experimental value of 1.17 eV. The final deviation 
is comparable to that of typical pseudopotential calcula- 
tions, which tend to a slight overestimation.^'^ Although 
a certain discrepancy with respect to the pseudopoten- 
tial results still remains, we conclude that it is, in fact, 
smaller than suggested by previous calculations. 
The stronger underestimation in these studies must be 
attributed at least in part to an incomplete convergence 
with respect to the number of bands (e.g., 24 bands in 
Ref. 9) in combination with the linearization error. In 
the case of Si we find that the latter accounts for less 
than 0.03 eV, but we cannot rule out that it is larger in 
other systems and must be taken into account in order 
to obtain reliable GW results. 

Very recently van Schilfgaardc et al?^ also reexamined 
the convergence of the self-energy and the quasiparticle 
energies in Si, employing the same GW algorithm^^ as 
in this work. However, the eigenfunctions were gener- 
ated by the full-potential LMTO method with basis sets 
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FIG. 4: Indirect band gap of Si as a function of the number 
of bands. The inchision of second-derivative local orbitals for 
I < 3 (B) provides the largest step towards convergence with 
respect to the basis set. 



ranging from 50 to 185 orbitals, including additional MT 
orbitals located in the conduction-band region. They ob- 
tained a T25'v-Xic gap of 1.15 eV, similar to our results 
in Table II; the slight discrepancy is due to the different 
way of generating the basis functions used to construct 
the self-energy. With the corrections for k-point conver- 
gence and screening due to the 2p electrons, their best 
estimate for the r25'v-Xic gap, if the GW approximation 
is evaluated with LDA eigenstates, is 1.10 eV. Further- 
more, Shishkin and Kresse reported quasiparticle band 
gaps obtained with a new PAW implementation in quan- 
titative agreement with our results. 



CONCLUSIONS 



band gap increases by less than 0.03 eV. This makes the 
LAPW method a suitable quantitative reference scheme 
for all-electron GW calculations. A stronger basis-set de- 
pendence in other systems or in the correction of higher 
bands is possible, however. If convergence with respect 
to the number of bands is also taken into account, the 
GW band gap turns out to be larger and closer to pseu- 
dopotential values than previously reported all-electron 
results, although a certain discrepancy still remains. 
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APPENDIX: SCALAR-RELATIVISTIC 
EQUATIONS 

The radial scalar-relativistic Dirac equations for the 
large and small components p and g of a free electron in 
a spherical potential at energy e are given by 



p'{e,r) = 2M{e,r)q{e,r) + -p{e,r) , 



(A.1) 



q'{e, r) = - -q{e, r) + w{e, r)p{e, r) (A.2) 

with M(e,r) = 1 + [e-V^ir)] / {2c^) and w{e,r) = 
+ 1)] / [2M(e, r)r2] + (r) - e, where y'=^(r) is the 
spherical part of the Kohn-Sham effective potential. 
Their i^th energy derivatives are given by 

pM'(e,r) = 2M(e,r)gM(e,0 + V'^He^O 

r 

+ 4g('^-i)(e,r), (A.3) 



In this work we investigated the influence of the basis- 
set accuracy in linearized all-electron methods on the 
GW self-energy correction. We showed that the addi- 
tion of local orbitals defined as second and higher en- 
ergy derivatives of solutions of the radial Schrodinger (or 
scalar-relativistic Dirac) equation constitutes an exten- 
sion of the LAPW basis set that allows a systematic im- 
provement towards basis-set completeness. In the case of 
silicon basis-set convergence of the GW results was essen- 
tially reached with the inclusion of second-derivative lo- 
cal orbitals. It is the basis-set dependence of (S'^) that is 
responsible for changes in the quasiparticle energies. As 
all of them decrease towards basis-set completeness, the 
effect on the relative transition energies is smaller than on 
the quasiparticle energies themselves. The fundamental 



gH'(e,r) = — q^""^ 
r 



(e, r) -I- ■u;(e, r)p^^^ (e, r) 



(i)(e,r)p('^-i)(e,r) 



(A.4) 



+ 



E 

H=2 



{e,r)p^-'-^\e,r 



where M'^^)(e,r) = l/(2c^) has been used. The uth en- 
ergy derivative w^"^ (e, r) has the expression 



w 



(e,r) = (-ir 



v\l(l + l) 



+ 1t 



Si. . (A.5) 



In the nonrelativistic limit these formulas correspond to 
Eq. (12). 
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